Roughness of Interfacial Crack Front: Correlated Percolation in the Damage Zone 
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We show that the roughness exponent C, of an in-plane crack front slowly propagating along a 
heterogeneous interface embeded in a elastic body, is in full agreement with a correlated percolation 
problem in a linear gradient. We obtain C, = v jiX + v) where v is the correlation length critical 
exponent. We develop an elastic brittle model based on both the 3D Green function in an elastic 
half-space and a discrete interface of brittle fibers and find numerically that v = 1.5. We conjecture 
it to be 3/2. This yields C, = 3/5. We also obtain by direct numerical simulations C, — 0.6 in 
excellent agreement with our prediction. This modelling is for the first time in close agreement with 
experimental observations. 

PACS numbers: 83.80.Ab, 62.20.Mk, 81.40.Np 



An important motivation for studying interfacial crack 
pinning |^] is to simplify the study of the origin of the 
scaling properties of brittle crack surfaces These 
scaling properties are seen, for example, in the height- 
height correlation of the fracture roughness (i.e. out- 
plane roughness), which shows self-affinity. That is, the 
conditional probability density /i), i.e. the probabil- 
ity that the crack surface passes within dh of the height 
h at position x when it had height zero at a; = 0, shows 
the invariance 



X''p{Xx, \''h) —p{x,h) , 



(1) 



where C is the roughness exponent. It is now generally 
believed that the roughness exponent shows a universal 
value of about 0.80 at larger scales while a lower value 
of about 0.5 might be seen at smaller scales [||. 

Direct observations of the interfacial crack-pinning 
have been performed recently. The problem consists of 
following the roughness of a crack front moving along the 
flat interface between two elastically connected blocks. 
The experimental study of constrained crack propaga- 
tion between two sintered plexiglass plates presented in 
resulted in the estimate of the in-plane roughness ex- 
ponent: C, — 0.55 ± 0.05. This work was followed by a 
longer study leading to the estimate ( = 0.63±0.03 0. A 
recent study of the motion of a helium-4 meniscus along 
a disordered substrate — a problem closely related to the 
motion of a crack line — gave C = 0.56 ± 0.03 B. 
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Numerous models for interfacial crack propogation in 
heterogeneous material have been proposed. The numeri- 
cal simulation presented in is based on a perturbative 
Green function approach following the quasistatic evo- 
lution of the interfacial crack front position a{x,t) [ p^ , 
which is treated as a single- valued function of position x 
along the orthogonal direction to the crack advancement 
direction, and time t. The linearized Green function used 
binds together points only along the fracture front. The 
stress intensity factor at a point x along the fracture front 
is then found to be 
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a{x' , t) — a{x, t) 
{x — .t')^ 



dx' , (2) 



where Kq is the stress intensity factor that would result 
if the crack were straight . The fracture is advanced 
by identifying the most stressed point along the fracture 
line and advancing this by a small step. The roughness 
exponent of the crack front was estimated numerically 
to be C = 0.35, while a direct dynamical renormaliza- 
tion group calculation gave C = 1/3 to lowest order 
Higher-order corrections to this quasistatic analysis in- 
creases the value of C to 0.48 iQ, while a different qua- 
sistatic analytical technique suggests 0.390 ± 0.002 [13|. 
Dynamic effects have beeen largely studied numerically 
and analytically jl^ , in particular in the form of crack 
front waves. They lead to an increased roughness expo- 
nent compared to the initial 1/3- value up to ^ = 0.5 but 
still smaller than experimental observations. 

Hence, the situation today is that there is a large and 
significant discrepancy between theoretical and experi- 
mental estimates: theoretical estimates being consider- 
ably smaller than the experimental ones. 

In this letter we present a numerical calculation of 
the roughness exponent of a crack front that propagates 
quasi-statically along a heterogeneous interface, based on 
a Green function technique which differs from that pre- 
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FIG. 1: The crack front for a 128 x 128 system. The fracture 
is propagating from bottom to top. The broken springs are 
black dots. The crack front is drawn as a white hne. 



viously used. Indeed it does not reduce the crack tip to 
a single tortuous line but describes the tip as a region 
of interactions between microcracks (see Fig. Our 
technique is based on the static solution of the elastic 
equations for the deformation of the surface of an elas- 
tic half space The local deformation, u, at position 
(x, y) along the plane is related to the normal stress field 
a by the expression 



G{x-x',y~y')a{x',y') dx' dy' , (3) 



where the Green function is llf 



Gix,y) 



(4) 



with s is the Poisson ratio and e the elastic constant. 
For large r = \{x — x',y — y')], G ~ 1/r. For the sake 
of simplicity, as generally done for the study of contact 
between two elastic bodies , we substitute one of the 
elastic plates by an infinitely rigid one. The other plate is 
modeled as an elastic block for which Eq. is valid. We 
discretize the model, so that the forces and deformations 
are described by the discrete version of Eq. O), 



(5) 



Gi.j is the Green function Eq. (Q) averaged over an area 
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+h/2 /.+6/2 



dx' dy'G(i-x',2-v') , (6) 



h/2 J~hj2 



where 6 is the lattice spacing, ui is the deformation of 
the elastic body at site «, and the force acting at that 
point. The indices i and j in Eq. ^ run over all I? sites. 

The elastic block is connected to the infinitely stiff 
plate by a discrete interface made of an array of elas- 
tic harmonic springs. Each spring is brittle and has a 
breaking threshold randomly drawn from a uniform dis- 
tribution between zero and one. The spacing between 
the springs is h in both the x and y directions. The force 
/,; that an unbroken spring i is carrying, is transferred 
over an area of size 6^ to the soft surface and given by 
Hooke's law: 




/i = -h{ui - D) , 



(7) 



<a>/L" 



FIG. 2: Family- Vicsek scahng of the crack front roughness: 
w/LP-^ vs. {a}/L°-^ for different system sizes: L = 16, 32, 64 
and 128. 



where k is the spring constant (k = 1 for all springs). D 
is the displacement of the infinitely stiff medium, and is 
a function of y, i.e. solid rotation of the stiff medium. 
The quantity (ui — D) is, therefore, the length that spring 
i is stretched. We assume periodic boundary conditions 
both in the x and y directions. In order to model the 
mode I fracturing of the interface between the two media 
(Double Cantilever Beam (DCB) test) in a way compati- 
ble with the biperiodic boundary conditions, we let D{y) 
be wedgeshaped with a maximum at ?; = and L, and 
zero for y = L/2 {i.e. mirroring conditions). The load 
the system, D{y) is increased and the springs break one 
by one. The numerical technique to solve the equations 
that ensue is presented in detail in ||l7|] . 

We show in Fig. |l| a typical damage front in a 128 x 128 
system. We define the fracture front a{x, t) in this model 
as the set of nodes that form a continuous path separating 
the infinite cluster of broken springs from the infinite 
cluster of unbroken springs. This definition is similar to 
that of Ref. j^. From Fig. |l|, we clearly see that even if 
a front can be defined, an extended damage zone exists. 
Accordingly the front does not capture all the active tip 
of the fracture. This observation supports the relevance 
of a non line description of the crack front. 

In F ig. PI we s how the width of the fracture front 
w = ^ (a'^J— (a)2 as a function of its average position 
(a) — which acts as a time parameter in this quasi-static 
model, for various system sizes. By collapsing the width 
evolution for the different system sizes, we see that the 
crack front follows a Family- Vicsek scaling |l^. Two 
important exponents are, thus, obtained: the roughness 
exponent C = 0.6 and dynamical exponent z — 0.9. The 
roughness exponent is in close agreement with the exper- 
imentally obtained value 0, |, D , while the dynamical 
exponent z was found in Ref. [|| to be 1.2. 

Our model distinguishes itself from earlier numerical 
models in three major ways: (1) Most previous models 
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FIG. 3: Fluctuations of the density of broken bonds Wc = 
((Pc) ~ iPc)^y^'^ plotted against L for a homogeneous parallel 
loading without any gradient, i.e. constant D over the whole 
system. The slope of the straight line is —l/v — —0.65. 



are based on a small perturbation approach and include 
linearizations pO[ | that are not used here. Furthermore, 
it was assumed in the earlier studies that (2) the frac- 
ture line was a single-valued function, hence ruling out 
overhangs. This assumption prevents islands of unbroken 
bonds from forming in the wake of the advancing fracture 
line. (3) The assumption of a single advancing fracture 
line also prevents the formation of a damage zone in front 
of the fracture line. None of these three assumptions are 
necessary in the present model. In order to test whether 
assumptions (2) or (3) are responsible for the difference 
in roughness exponent, we imposed both on the present 
model. No change in C was observed. Hence, we conclude 
that it is the linearization assumption that is responsable 
for the difference. 

In the following, we demonstrate that the roughness 
exponent C is a result of a correlated percolation pro- 
cess in a gradient imposed by the mode I loading of the 
system. We start by charactarizing the correlated per- 
colation. For this goal, we consider a similar problem 
but where the loading is obtained without any gradi- 
ent, i.e. horizontal and parallel displacement of the rigid 
medium [ pT[ . In this case, when the homogeneous dis- 
placement D of the rigid medium reaches a maximum 
value, Dc the system goes unstable, and unless D is de- 
creased again, catastrophic failure sets in. In Ref. 
the size distribution of clusters of broken springs was 
studied, and a power law was found with an exponent 
T = 1.6. This value is different from ordinary percolation 
where r « 2.05 |2^ and shows that correlated percolation 
takes place. In Fig. we show the fluctuations of density 
of broken springs a,t D = Dc. If there is a diverging cor- 
relation length in the problem, these fluctuations scale 
as L-i/"^. We find that = 0.65, leading to ly ^ 1.54. 
We conjecture that the exact value is v = 3/2. Hence, 
the fracture process in this system is in a universality 



FIG. 4: Damage profile p{y/L) in reduced variable y/L for 
three different system sizes L. 



class which is different from standard percolation where 
v = 4/3. 

When D is no longer constant, but is given by the 
wedge shape (DCB load), the form of interactions be- 
tween the springs does not change. Hence, the critical 
properties of the model when run in the "horizontal" 
mode are still present under DCB loading. The wedge 
shape of D leads to a gradient in the loading of the 
springs, going from very high loading where the dam- 
age is large to very low loading well into the still intact 
part of the interface. In Fig. ^, we show the damage pro- 
file p{y), i.e. the density of broken springs averaged in 
the X direction, for systems of different sizes. The profile 
is clearly linear. Somewhere along this damage profile, 
there is a line in the x-direction at y = where the dam- 
age density is critical, p{yc) = Pc- In the vicinity of this 
line, there is a critical region which is characterized by 
being on the edge between stability and instability and 
corresponds to the crack front. Following the arguments 
of Sapoval et al. for percolation in a gradient pT[ , if p{y) 
follows a linear law, 



p{y) = 1 



y-yo 
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(8) 



on the interval yo < y < yo + ly, where ly is the length 
scale characterizing the width of the damage zone, then 
yc = ya + lyi^ — Pc) is the position of the fracture front 
(when ignoring small corrections for finite-size systems). 
The critical region has a width ^ = \yw — yc\ which is 
related to the damage p{y) as ^ ~ \p{yw) — Pc]^" . Ehmi- 
nating \yc — yw\ between these two expressions for ^ gives 



Ly 



(9) 



Since the gradient in D is inversely proportional to L 
and the average strength of the springs does not change 
with L, we have that 



ly (X L 



(10) 
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In Fig. |j, the reduced variable y/L was used resulting in 
data collapse for different system sizes, thus validating 
Eq. (|ic|). Furthermore, the width of the fracture front, w 
is proportional to the width of the critical region ^ oc w. 
Hence, we find 

vjr^L^^L^/^. (11) 

where we have used our estimate v = 3/2. This result 
is in excellent agreement with our numerical simulations 
and with the experimental results. 

A similar idea has been proposed for the outplane 
roughness of fracture surfaces where the gradient and 
the V exponent are different and leads to C = 4/5 in ex- 
cellent agreement with the experimental value 0.80 [^ . 
We finally note that Zapperi et al. l2|] have pointed 



out the connection between gradient percolation and the 
interface fracture problem in the limit of infinitely stiff 
plates. We have shown in this letter that this mechanism 
stays intact also when the plates do respond elastically. 
However, the universality class of the corresponding cor- 
related percolation problem is different from ordinary 
percolation, leading to the observed roughness exponent 
C = 3/5 rather than 4/7 for ordinary percolation in a 
gradient. 
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